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ABSTRACT 


A Liquid Metal Fast Breeder Reactor is simulated on a hybrid 
camputer to study the transient behavior of the neutron density during 
a controlled reactivity input disturbance. 

The nonlinear partial differential equations of heat flow are 
reduced by a discrete space-continuous time method to ordinary non- 
linear differential equations which are readily solved on the analog 
computer. Use is made of time multiplexing of the analog circuitry 
in order to reduce the number of camponents. An open-loop iteration 
process is employed to solve the closed-loop feedback system. 

Recently conducted research with the open-loop iteration method 
has demonstrated that a large number of iterations are required for 
convergence. An algorithm is developed which gives an improvement in 
the convergence rate for early values of time. For times greater than 
two seconds a stability problem with the converged solution was en- 
countered and is discussed with some observations and comments. 

Innovations are introduced in the simulation of the neutron kinetics 
equations and in the handling of the nonlinear thermal properties of the 


uranium oxide fuel. 
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I. INTRODUCTION 


For over twenty years the nuclear community in the U. S. has been 
working to meet the recognized important objectives of economic nuclear 
power and improved fuel utilization through the development of the 
Liquid Metal Fast Breeder Reactor [Ref. 1]. This effort, however, 
has been largely dispersed with no real organized effort. Starting 
in 1965 an unprecedented demand for light water reactor power plants 
occurred paralleled with increased uranium needs. This increased fuel 
requirement emphasized the need to develop a Liquid Metal Fast Breeder 
Reactor and has caused the U. S. as well as other countries to establish 
Liquid Metal Fast Breeder Reactor (LMFBR) programs. 

The light water reactors are very successful today with favorable 
prediction of continued improvements in fuel performance and reduced 
operating cost [Ref. 2]. However, the LMFBR can provide the most ef- 
ficient means of exploiting the energy available in uranium because of 
a high thermodynamic efficiency, an efficient fuel conversion cycle, 
reasonable capital and operating costs, and safety in operation. 

The present U. S. effort in Fast Reactors is summed up in the 
nuclear reactors Clementine (0.025 Megawatts), LAMPRE I (1 Megawatt), 
EBR-I (1.4 Megawatts), EBR-II (62.5 Megawatts) and Enrico Fermi (200 
Megawatts) of which only the latter two remain in commission. The 
Current LMFBR program envisions a 1,000 Megawatt electric reactor, 
which is a substantial increase over the largest system constructed 
to date. It is generally accepted [Ref. 3 and 4] that the operating 
condition of the proposed LMFBR is beyond the range which permit prudent 


extrapolation from existing experiences. Simulation studies can greatly 


alba, 


asSist in design evaluation by predicting the behavior of the camposite 
model based on experimental information obtained from the component 
Dakus: 

The use of hybrid simulation is a very effective method to carry 
out design evaluation based upon the transient response to predictable 
accidental disturbances. The pure analog approach is not practical due 
to the prohibited amount of hardware required to simulate the coupled 
nonlinear differential equations of heat flow. The pure digital ap- 
proach, which is feasible, does not readily lend itself to multiple 
parameter studies for design optimization; the computer time required 
would be economically unfeasible. The hybrid computer combines same 
of the best characteristics of the analog and digital machines result- 
ing in a more economical computer. Recently, large-scale simulations 
[Ref. 5 and 6] were performed with the object of camparing the economics 
of hybrid versus digital for large plant simulation. These simulations 
indicate that this advantage is not about to disappear. 

Previous efforts have been made to demonstrate the feasibility of 
employing a hybrid computer for the analysis of the transient behavior 
of a LMFBR core [Ref. 7]. The development of new techniques to enhance 
the effectiveness of the hybrid computer are currently in progress 
[Ref. 8]. The method developed in this report is a contribution in 
this effort with regards to software design related to reducing the 
number of iterations for convergence. The number of iterations is 
Critical in a hybrid simulation since a large number of iterations 
would defeat the purpose of using the hybrid computer. A pure digital 
approach would then become more practical where both discrete time and 


discrete space could be employed. 
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Additional contributions were realized by simulating the nonlinear 
thermal conductivity of UO., on a Diode Function Generator resulting in 
a reduction of analog hardware, and simplification of the neutron 
kinetics equations for use on the analog computer. 

The parameters for the model used are primarily from a recently 
proposed 1000 MWe LMFBR design [Ref. 9]. However the model developed 
is quite general and can be employed in the analysis and design studies 


of fast, intermediate or thermal reactors. 


iy 


II. DESCRIPTION OF MODEL 


A. REACTOR KINETICS 

The point reactor kinetics equations are used to describe the 
neutron density behavior. The assumption is made that the density 
(F) can be separated into a time dependent and a space dependent 


term, and can be expressed as” 


pe Y(R) ne) (1) 
Reference 10 and 11 present good treatments of the separation of 


time and space of the neutron density. n(t) is given by the fol- 


lowing equations 
N 
By l4) 7 | pf) - Ely.) + Z » ce CE) (a) 


‘ (2) 
ee it (5, 
fo. Poi ee) (4) 
it A 
where 6-2 3 and with initial conditions /((¢)- re A Cen 
AL 


FC. 


A treatment of the point reactor kinetics can also be found in Ref. 
IO and. dee 

In most applications of the above equations it is considered 
appropriate to treat N = 6, that is 6 delayed neutron precursors. 
However, in this model due to lack of sufficient analog hardware only 
one delayed precursor will be considered; its parameters will be 


appropriate averages of those of the six. 


tR refers to the three spacial variables of the core and not to 
be confused with r, the radial variable of one fuel pin. 
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Equation (2) is recast into a more convenient form by defining 


two dimensionless quantities (see Appendix A for details) 


eoeyult) meercailt) 
Ia iad ares ae and Dt) = mes 


With these definitions and the one precursor assumption equation (2) 


becames 
ANA) P/V ety 
D(t) 
= = ” | Not}- D(t) | ‘ta 


with initial conditions N(0) = D(0) = 1, and where the term K(t) 
equals P (t ps and is defined in dollar units. p(t) is the sum of 
the forcing function P (t), plus the feedback reactivity is (t), 


or simply 


Pie = Rae) Pe) (4) 
The quantity P/ yi for a fast reactor is of the order Tor thus equation 
(3a) requires a high gain integrator for similation on the analog computer. 
This introduces some problems in the simulation which will be discussed 
in Section V. If equation (3a) is rearranged as follows, 
KNEE: NE) + D442) - (YE) NG) 
N (2) represents a small quantity campared to the other 


bed 1) 


O 
the term (44 
three terms of the equation. Neglecting this term is a reasonable ap- 


proximation and equations (3a) and (3b) become 


N(th= K(t) NG) + BCt} (a) 
(5) 
Be = |] m(t)-ddt)) (b) 


1 


Equation (5b) is an initial value ordinary differential equation and 
its solution is eaSily handled by the analog computer. 

The forcing function is an arbitrary selected but known function 
of time. In this report a terminated ramp is used to model the with- 
drawal of a control rod. The feedback reactivity is not a known 
function of time but depends on the average temperature of the reactor. 
In a fast reactor there are three main contributions to the feedback 
[Ref. 12], the Doppler effect, fuel expansion and coolant expansion. 
The expression for feedback due to these three effects can be repre- 


sented ag 


ital 


haj= ms soon, An(\+ als rd Byes les map orm AT. 


1s 


or, in terms of dollar units 


z eee. pea — 


E é 


eygeh a and C Kool are positive constants. The temperatures 
indicated in equation (6) represent average core temperatures. 

In this model a single average fuel cell is assumed to give all 
the temperature information necessary to construct the average core 
temperature.” 

The power per unit volume supplied by the fission process can be 
expressed as 


Q = Nit) SCR) 


where S(R) is proportional to \y(R) . The proportionality factors include 


oe Appendix A for more details. 


References 8, 13 and 14 all employ a single "average" fuel cell 
to compute the reactivity feedback. Reference 7 uses five fuel cells 
properly weighted to compute the reactivity feedback. 
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eO\p E (average fission cross section), E (energy released per 
fission), and V (the average velocity of the neutrons). In this 


model E j > ¢ and V are treated as constants. 


B. HEAT FLOW 

Large power densities exist within the core of a 1000 MWe LMFBR 
giving rise to steep temperature gradients and extreme heat flow. 
Furthermore the thermal conductivity of the ceramic fuels being 
considered for the 1000 MWe LMFBR (namely, oxides, nitrites and 
carbides) vary nonlinearly at the proposed high operating temperatures. 
These factors dictate that the problem of heat flow be approached on 
a microscopic basis. 

A LMFBR is composed of thousands of pencil size fuel pins approxi- 
mately three feet long in the shape of right cylinders; the coolant 
flows coaxially. A sample fuel pin illustrated in Fig. 1 is composed 
of a cylindrical shaped fuel assembly with a concentric shell called 
cladding. 

The heat flow equation for a cylindrical shaped body with varying 


thermal properties is: 
A Looe B] 6 2 bene 


o [cong || = QU, 2,2,t) 


In a fast reactor the neutron flux is nearly constant within the fuel 


ey 
C = 
(8) 


pin in the radial direction (in contrast to a thermal reactor where flux 
depression occurs within the fuel). This fact gives rise to cylindrical 
symmetry and if axial flow of heat is neglected, which is consistant 
with other models proposed [Ref. 5, 13 and 14], equation (8) for any 


fuel pin reduces to: 


2a 





rn 





—}—-—} ae TI 


CO OO Oe ET TT A ET TE EE Sy = ao 


aa et NY TOO > 


: 


a ee eee ee ee —— a ee ee ee ee SE STEERER To Eee ft 


NIGG TO 








Diagram of Fuel Pin 


Figure 1. 
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Fuel:  ¢, 4; ae + es | v K, (7) | 4 $(2) N(t) (9) 


dv Ov 
| iret o ve 
Cladding: (¢, d, eae | y Ke (T) ai (10) 


where the assumptions outlined in Section IIA for the power generation 
Q, have been introduced (i.e, the space and time separation of neutron 
density). The term S (R) has been replaced by S(z) to indicate the 
heat generation rate in a single fuel pin is a function of the axial 
direction only. 

Since the analog computer can handle only ordinary differential 
equations, equations (9) and (10) are reduced to initial value ordinary 
differential equaticns by employing a method of integrating over the 
spacial variable [Ref. 15]. The technique involves discretizing the 
space variable in the radial direction, Y, into five nodes within the 
fuel and one within the cladding (see Fig. 2). In carrying out the 
integration the temperature is assumed to vary linearly between nodes. 


For more details of the derivations refer to Appendix A. 


As a result of employing this technique equations (9) and (10) became: 


Fuel? €.d A (3 Py | 7 bg 
are (s 24V; ae 74 Ts (3 He au pre 
K ae a! _ T a 
7 | oY Lt) = 
at Fiitiie (|+ ZX. \( are - 





STVe ie) (i- = am aa is) s(2) NQ) 


(1 indicates the node which varies in number from 1-5 within the fuel) 


Kit — a conductivity at boundary between node i and 
i 
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— or ee + eee ee ee 


It~ 


Kradiuc 


Temperature Profile Assumed Within Fuel Pin. 


Fagunces. 
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Cladding: Cede a AS. (q- Ss )+ Sime (d+ =.) 


“5 (12) 
D sig tert ys ; 
axe /2. 


ae. 
=e Je 
Ke ( | 2¥.) ( x | 
(Superscript j refers to the axial node) 


Due to symmetry about r = 0, the condition that We = a is applied 


as a boundary condition. Also at the clad-fuel interface the continuity 
of heat flow requires that 


aie I 
Kia E : ‘i K.¢ 


Fas VA 








At the boundary between cladding and coolant the conservation of 
energy is applied to arrive at the following equation (see Appendix A 


for derivation). 





dk , G4) gk _ _24 
ot ae 2 7 we a Bi cae (13) 


G(t)= mass flow rate of coolant 
H = heat flux across the clad-coolant interface 


r_ = radius of fuel pin (fuel + cladding) 


Again an ordinary differential equation is needed. Here, however, 
the use of finite differencing of the space variable is more practical. 
This simply requires replacing Teh by Ole ll Wh S2 and assuming 
the coolant temperature varies linearly over length AZ. Thus equation 


(13) can be expressed as: 


ATS Ct) (T: a ie - Za 
. ~ ¥ 





Ht dB 


a ' es 5 
"a =average coolant temperature of a axial segment = <= i 
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At the clad-coolant boundary the heat flow per unit time (H) 
is represented by h(Ty- T. ) which also equals K, ( and, : 
h is the heat transfer coefficient between clad and coolant. 

After incorporating the boundary conditions and the physical 
constants listed in Appendix C the following coupled ordinary dif- 


ferential equations are obtained. 


1 -_ 3 _§ — 
ae [T'+, 26337; | = 328.2 ]K 1. (Te J+ .4) $42) N (t) 


(15) 
dos = -—7 +I 
a | Te) +, 203m omegs: Th! |e 23/.0[7K,,,(G-h) | Ts 
ae (TT) ip, 433 $G@) N() 
. Tir 1297 + A44T Je 208.0 [K,430G-B 
‘ i7 
. /3BWeMK, 5; ( Tylor |) 2 933 52) Ai Ct) (17) 
ae ¥ _ ner 
Tt | T+, 1826 lot. 1507 ty; [= 198.1 ee. (Ts'-T) | (/¥) 
- 1485 1K, 43 (T-T)] +. 433 82) NU) 
i | —) a? 
At as + ./840T% +4362 Io ]e 464.0 LK ¢ (1'-75') | 
ae } ; 
se Ear ah -T? J + 528 s02)NQ) a 
ji | | =n 
ie [T+ 1.3257) +.aesT?! Je 174.1 (G-T!] 
za 
mG? .o L Kye Co Teal ii 
Ae a a 
G, let — 166.37! + 18.8% + 6757, (21) 
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The inlet temperature to the first axial segment is known and 
treated as a constant during any transient. In nuclear reactors 
the inlet temperatures are usually kept constant, and any rapid 
power change will have leveled off before any change occurs in the 
inlet temperature. For each subsequent segments the inlet temper- 
ature is a known function of time obtained from analysis of the 
previous section (see Digital Software section for details). The 
Steady state spacial power distribution S(Z) is known a priori, 
therefore, equations (15) through (21) can be solved on the analog 
Sempucer . 

The bonding between the fuel and cladding is ignored in this 
model and the fuel is assumed to be in constant contact with the 
Cladding. In most LMFBR under design the bonding agent is helium. 
Neglecting this bonding will produce lower than normal temperatures 


within the fuel. 


Jag 


III. HYBRID COMPUTER TECHNIQUE 


A. DESCRIPTION OF ANALOG COMPUTER HARDWARE 

From the heat flow equations given in the previous section it can 
be seen that the model requires the simulation of a large number of 
ordinary differential equations. However, an important aspect to keep 
in mind is that the description of each axial segment incorporates 
the same differential equations, the coolant inlet temperature and 
the fuel heat generation rate being the only quantities varying fran 
segment to segment. With this in mind a technique of multiplexing 
that is, time sharing, of the analog circuitry can be employed. The 
details of the multiplexing used will be explained in Section IIIB. 
To use the technique, only one set of equations (15) through (21) 
need be set up on the analog computer. In order to solve equations 
(15) through (21) on the analog computer they must first be magnitude 
scaled to insure that the outputs do not exceed the maximum analog 
computer voltage (+ 100”). After proper scaling equations (15) 
through (21) can be expressed as equations (41) through (47) which 
are shown in Appendix B. 

The approach used in setting up the solution to equations (41) 
through (47) on the analog computer can best be explained by taking 


one equation, 


rt |T +.2038e T+. nas! ]= 23.1 fr Ki 2,3 T3-7)] 


ee Se (42) 
— 1.55 [10 Ky 12 Cy -T!)] + 2165 s@ NCE) 
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The right side of equations (42) represents inputs to an 
integrator whose analog symbol is shown by element A in Fig. 3. 
Element B is a summer amplifier. This process is extended over 
the remaining heat equations. Figure 4a is a diagram of the analog 
Gryecuirery. 


SS 


LCs te + 20387 +. 12957, 


29528 Ts. 








ee 


=. 716 7, +.2038 Ty" Pees 7” 





~ ae 


I.C. refers to initial conditions. 


Figure 3. Illustration of Technique Used to Construct Analog 
eared t . 

The nonlinearity of the thermal conductivity, K(T), was readily 
handled by the analog camputer through utilization of a Diode Function 
Generator (DFG). A DFG can be used to approximate a nonlinear curve 
by the method of straight line segments. The DFG used with this 
simulation has ten break points. Thus a maximum of ten straight 
lines were available to duplicate the curve and this was quite adequate 
to give a good representation of the nonlinearity. Experimental data 
for the K(T) curve was obtained from Ref. 16 and 17. Figure 5 shows 
both curves. The DFG operates simply as shown in Fig. 6. The input 
is a voltage which represents ok temperature and the output is a volt- 


age proportional to the thermal conductivity for that temperature 
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Figure 4b. Diagram of Multiplexed DFG. 
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Figure 4c. Diagram of Logic Circuit. 
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Figure 4d. Key For Analog Symbols. 
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Figure 5. Thermal Conductivgty Of UO.» 
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Figure 6. Block diagram of DFG. 


The boundary conditions of the five fuel nodes required five 
values of K(T). However only four DFG's were available on the analog 
computer, consequently, one of the DFG's had to be time shared. The 
time sharing was accomplished by alternating the information to a 
DFG fram two boundaries. The switching was accomplished by two 
digital/analog switches. The logic to operate one of the switches 
was supplied by a 10 kHz square wave which for clarity will be called 


logic variable L,; the other switch was operated by Ly the inversion 


il 
of Ly. The output of the DFG was fed to two track-hold amplifiers, one 


being controlled by logic L. and the other by L,. With the aid of Fig. 


a ah 

4b it can be seen that while one switch is closed one amplifier will 
be in the track mode following the output of the DFG while the other 
is in the hold mode. When the logic signal inverts the switches 
reverse conditions and the amplifiers reverse modes. 

The analog computer used for this model, COMCOR-5000, is equipped 
with digital logic and in addition completely controllable by the 
general purpose digital computer, SDS-9300. However it was necessary 
to develop a logic network which could control the compute mode of the 
analog computer automatically using the digital computer for varying 
periods of time. Furthermore at times it was necessary to bring the 
model up to operating temperature before the external forcing function 


was applied to the reactor kinetics. To realize this flexibility a 


complete logic network was developed which is shown in Fig. 4c. 
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B. DESCRIPTION OF DIGITAL COMPUTER SOFTWARE 

In the previous section it was pointed out that the solution to 
the heat flow for the complete fuel pin involved one set of equations 
by applying the analog circuitry of Fig. 4a on a time sharing basis. 
The digital computer with its logic capabilities can control this time 
sharing very effectively. Essentially what occurs is that the analog 
computer is used as a large subroutine and when called produces the 
solution for the axial segment under consideration. 

In this simulation the solution to one axial segment is treated 
on a continuous time basis during any transient condition. In order 
to treat each segment in this manner the time behavior of the inlet 
coolant temperature must be incorporated within the solution for that 
segment. The digital computer with its memory capabilities can ac- 
complish this effectively. The digital computer through an analog- 
digital converter samples at a fast rate (250 Hz) the outlet tempera- 
tures of the segment being considered. This information is stored as 
an array within the digital computer and used as input data in the 
form of a analog voltage through a first order digital to analog 
converter during the solution of the next segment. Thus the solution 
can continue sequentially from one axial segment to the next. 

The feedback shown in equation (6) is an algebraic function and 
for this reason it iS more convenient to have the digital computer 
handle this aspect of the simulation. However, it is necessary at 
this point to note that the feedback mechanism is a function of the 
average temperature. This situation requires that the average tem- 
perature be known at a given time in order that the digital computer 


can evaluate the feedback. However the average temperature of the 
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core is not known until all the segments have been considered. To 
obtain the average temperature, the average fuel and coolant ten- 
peratures of each segment are weighted to approximate the average 
temperature. Each segment is then evaluated on a continuous time 
basis with the digital camputer supplying the closed-loop feedback. 
The operation is illustrated in Fig. 7 where the closed-loop system 
for each axial segment is shown. During the analysis of each segment 
the digital computer rapidly samples the logic condition, defined L3, 
on a 250 Hz signal which commences simultaneously with the signal 
which places the integrators in the campute mode (Fig. 8). When L3 
is -6 volts the digital computer samples the change in ambient tem- 
perature for the fuel (AP?) and coolant (AF?) for the axial segment 

} being considered, stores the information, camputes a feedback 
after weighting the temperatures, and outputs the feedback just 
computed as a voltage to the analog computer (Fig. 8c). L3 during 
the above operation has returned to 0 volts. The digital computer 
remains in a loop until L3 again has the value -6 volts whereupon 
the process in Fig. 8c is repeated. The condition of O volts on 
logic signal L2 occurring at time t ends the camputation for that 
Segment. For clarity the analysis of all ten segments in this manner 
will be referred to as iteration l. 

After each segment is analyzed the fuel and coolant temperatures, 
each represented by a 250 x t matrix, are stored on a disc due to 
storage limitations within the digital computer. After all ten seg- 
ments have been analyzed the 10 x 250 x t array for each temperature 
1s processed to obtain a 250 x t array representing the average tem- 


perature as a function of time. A K, (t) is then computed based on 
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this average temperature. A new iteration is made using the K, (t) 
just computed. The changes in fuel and coolant temperatures are 
sampled as before. However this iteration differs from iteration l 
in that the loop is not closed but open. Figure 9 shows the sche- 
matic for the open-loop process. At the conclusion of analyzing the 
ten axial segments a new K, (t) is computed based on the temperature 
profiles just developed. The open loop process is then repeated 
until convergence of the neutron density N(t) occurs. 

When applying the open-loop iterative method, a study should be 
conducted to insure that convergence will occur before attempting 
this procedure. Proof that this system will converge is given in 
Ref. 5 and 18, however Section V of this report covers a discussion 
on the rate of convergence. 

Figure 10 represents a flow chart for the digital software program 
just outlined. 

The method for obtaining the weighting factors for the coolant and 
fuel temperatures is now described. Normally to start the problem 
from a known steady-state spacial power distribution, the initial con- 
ditions on the itegrators were needed. These were easily obtained by 
putting in a step input of heat and allowing each segment to came to 
Steady state (operating level). The final voltages on all integrators 
were sampled and the information put on paper tape for future use as 
initial conditions. After each segment had reached operating level, 
it was disturbed with Kt) as a ramp input for a short period of time 

.l sec) without feedback. The changes in average fuel and coolant 
temperatures for this time interval were sampled for each hk axlal 


segment and from these temperatures the changes in average coolant 
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(aT) and fuel (AT) temperatures for the fuel pin were calculated. 


The weighting factors, Wy and Ww, were calculated using the following 


relationships, 
| a: 
ame 
Wie Sk 
oat, 
These weighting factors were then treated as constants for the 


duration of the transient condition. 
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IV. DYNAMIC CHECK OF MODEL 


When a physical system is described by a model it is important 
to know whether this model is adequately describing the system. 
To ascertain this, parts of the mathematical model were compared 
with known mathematical solutions or with digital computer solutions. 
The approximated neutron kinetics equations (5) were compared 
with the exact solution for the case of the one delayed emitter model 
using a step input of reactivity of magnitude 50 cents. Representing 
this step input by Ky the solution of the reactor kinetics can readily 
be obtained for this input condition [Ref. 19] and is given in equation 


(22) for the one delayed emitter model. 
Bae at ce bt 
Nt) = ( meer iu ( oy K (22) 


where, K » 


When practical numbers are used in equation (22) product terms con- 


taining A can be neglected and the solution can be expressed as 


} ( 
Nit)= izK 2 "a, (22a) 


Since 4 re ioe the second exponent has a time constant 2 x 10% 


sec., therefore it is unnoticable on any time scale of interest and to 
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a good approximation it can be neglected. What occurs with this 
assumption is that N(t) is a step jump to "a instead of 
delaying ~ .001 sec. Figure 11 shows both the analog solution of 
equation (5) and the exact solution given in equation (22a). From 
the results it can be seen that the approximation used to develop 
equation (5) appears valid. 

Equation (5) was also solved on the analog computer using a 


triangular pulse input for K(t) as follows: 


-5t dollars, of t © | sec. 
K (t+) = (i1-.5t) dollavs, lsec. <t © £ SRC. (22b) 
O Pee 25S G: 


) 


A digital camputer integration solution utilizing the Runga-Kutta 
method with a Adams-Moulton predictor corrector with error check 
(absolute error <¢ 107 per integration step) (Ref. 20 and 21] was 
also used to solve the reactor kinetics without the approximation, 
equation (3). The results of both solutions are shown in Fig. 12. 
Table I shows a comparison between the values for the two solutions 
for every tenth of a second. The agreement is quite good. 

The heat flow equations were solved on the analog computer for 
a homogeneous material with constant thermal properties. The exterior 
of the material was kept at a constant temperature and a constant heat 
source was applied within the material. The solution for this situ- 


ation, recalling that the equations apply for a circular cylinder, is 





[Ref. 22] 
K 2 
qi@-r') 2g 2 a i, oa 
T(t) = 4K aK Wer 2 xX, J, (K&4) 
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Figure 11. Camparison of Analog And Exact Solution For Neutron 
Kinetics For a Step Input of 50 Cents. 
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Figure 12. Camparison of Analog and Digital Solution of Neutron 
Kinetics for a Triangular Pulse Input for K(t). 
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TABLE I 


COMPARISON OF DIGITAL AND ANALOG SOLUTION OF REACTOR KINETICS 


Time Digital solution Analog solution 
Ora 120154 12052 
0.2 1. ee J, Ey. 
O23 1. ia 12192 
0.4 Ase L282 
O25 1293 1.390 
0.6 e525 i. 322 
On 1.689 1.685 
0.8 1.895 Loon 
0.9 2.157 2 a0) 
alo 2 S00 2a4eu 
i ngewe 2.398 26306 
Lee 2e205 20205 
re: 2ao4 2 
alee 2.095 2.075 
ces 1.999 1.980 
IPSC 1.906 1.890 
7 sleulis 1.800 
Wes 1730 Lae 
1.9 1.646 1.629 
2.0 1.566 12556 
Pll OCG PB EEK 
20 1.566 1.558 
uae 1566 2556 
ee 2566 13558 
Des 1.566 12556 
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where, 


heat applied in Watts/em™ 


C= 
a = radius of cylinder 

K = thermal conductivity of material 

eae ee root of Bessel function of order zero. 


Figure 13 shows a comparison between both solutions. It can be seen 
that the solution is quite good for nodes 1 through 4 but node 5 is 
in error by as much as 25% in early values of time and 10% in final 
time. The situation is generally recognized [Ref. 23] as being a 
problem near abrupt boundary conditions. Since abrupt changes are 
not envisioned as part of any study with this model it can reasonably 
be expected this error to be less. Furthermore the average temper- 
ature is the important quantity of interest; thus the error will be 
further reduced in magnitude. 

Prior to making any computations with this model it is wise to 
carry out a static and dynamic check of the model to insure all 


components are functioning properly. 
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V. RESULTS 


The principle objective of this report was to simulate the LMFBR 
on the hybrid computer and to improve on the rate of convergence in 
the open-loop iterative method. In the development of the technique 
used for improving the rate of convergence it was discovered that a 
problem related to the stability of the converged solution existed. 
Because of this problem of stability of convergence (to be covered 
later in this section) the solution of the problem is limited to a 
maximum of 2 seconds. 

Figure 14 and 15 show the N(t) response as a result of an 
external disturbance of one dollar per second which is terminated 
at the end of one second. The disturbance is intended to simulate 
the withdrawal of a control rod. The disturbance was started fram 


a steady state operating power level of 


sq) = (1000 Watts /cw? ) sim (Tm 

L represents the length of the fuel rod which is 100 cm; z is the 
axial position variable. Figure 14 illustrates the process of 
convergence of the open-loop iterative technique. Iteration number 
l was a result of a guess for the feedback. Ten iterations were 
required to obtain the converged solution. Figure 15 shows the 
response of N(t) produced after the iteration which uses the weighting 
factors, and this compares favorably with the converged solution 
shown on the same graph. 

During the development of this model it became questionable 


whether the open-loop iteration scheme would produce a converged 
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solution. Even though there exists mathematical proof that con- 
vergence occurs [Ref. 7 and 18] and published results [Ref. 7] 
showed the method successful, early results with this model re- 
vealed a divergent tendency for values of time greater than 2 
seconds. 

In an effort to isolate the problem a simple integrator circuit 
(Fig. 16a) was considered, since it can be shown mathematically that 
convergence is assured when using the open-loop iterative process 
[Ref. 18]. For the open-loop iterative process, the circuit in 
Fig. 16a can be schematically represented by that in Fig. léb. The 
procedure followed was to study the behavior of the open-loop iter- 


i representing the closed-loop analog feedback 


(2) 


ative process with Y 
for the first iteration. The output Y was sampled by the digital 
camputer and used as input for the succeeding iteration. This process 
was continued and the following observations made: 

a) Starting with a "good solution" for y (2) the behavior of E, 
with successive iterations was as shown in Fig. 17 (greatly exaggerated). 


(2) 


The region where Y began to diverge from the converged solution was 


when the output approached steady state condition, or when the sum 
On cs ym. The output would eventually converge but for times 
greater than 4 seconds (time constant for this example was 1 second) 
it tended more to oscillate with approximately equal error above and 
below the converged solution. 

b) For times less than four seconds a converged solution occurred 
but it could not be maintained indefinitely. Periodically (with no 
apparent consistency) the converged solution would be lost. Figure 18 


shows a converged solution and a succeeding one to illustrate this point. 
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Figure 16. Diagram of Integrator with Feedback. 
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The system would undergo a convergence process to again obtain a 
solution. Periodically (again with no apparent consistency) while 
in the process of converging it would lose the tendency to converge 
fram one iteration to the next and the process would again have to 
repeat itself. For times up to about 3 seconds, when X + y (2) was 
large, the output appeared unaffected and maintained a converged 
solution for repeated iterations. 

c) The rate of convergence was slow for interval of times greater 
than 4 seconds. 

To study the behavior of the neutron kinetics equations using 
the open-loop method a Simple lumped-parameter model of a nuclear 


reactor was developed. The mathematical model for the reactor is 


given by equation (25). 


Nit) = KG) N(t) + DY) (4) 
pit) = ALNU) — DUD) (b) 
=i = aT + ANC) (e) Wey, 
K(t) = Ke (t) + K;, (4) (d) 


_ Avese [AT 
Kt) = - - aes (e) 


where, T = average temperature, 
AT = change in average temperature above Ty the initial temperature, 
1/a = time constant of fuel. When heat is removed the temperature 
of fuel will fall to l/e of its steady state value in time 
equal to the time constant, 
A = factor converting neutron density to power. 
The Doppler effect was the only feedback mechanism considered and it was 
approximated by equation (25e); a closed-loop analog solution could then 


be obtained. Again the above equations can be schematically represented 


2 


as Fig. 16b. The open-loop iteration solution was attempted using 
(1) 


the analog solution for Y Similar behavior to the observations 
made above was noticed. 

Two factors to be considered when using the open-loop method 
based on the above observations are the rate of convergence and the 
Stability of the converged solution. The rate of convergence for the 
system described in Fig. 16 is slow. This can be shown as follows. 
If X (Fig. 16b) is considered as a step input of magnitude one, then 
for an initial guess of yh) (t) = 0 the output, y (2) (t), is equal to 


Gt (recalling the gain of the integrator is G). For a unit feedback 
ee 


the input for iteration 2 is 1 - Gt, and the output is Gt - oe, Tab 
is easy to show that a continued iteration process will develop the 
following series, ; P “ | 
z 3 4 L+/ Eye 
vite ce- Gh + Chae oe. ei 
= Z , u! 
where i denotes the iteration. A necessary condition for the series 
re 
to converge is lim es © - As long as G and t are finite this con- 
jroo ¢: 
dition is met. However before convergence begins to occur the ratio 


of the i+ 1 term to the jth term must be less than one, or 


— <i (26) 
The value of G used in Fig. 16 was 10. Thus convergence can be seen 
as a Slow process. It is also possible to understand the initial 
divergent tendency of the open-loop method conmented on earlier. 


Recalling that y 


in Fig. 17 represents the closed-loop system 
solution used for the first iteration, then if their exists any dif- 
ference between the step input X of Fig. 16b and the feedback y 
near steady state the integrator will continue to show a change in the 


output with the later values of time associated with the larger 
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differences. If it is assumed that any difference between the analog 
solution and the digital solution is constant in time, it can be seen 
fran the criteria used in developing equation (26) that the solution 
for early times will begin to converge while later values will tend 
to diverge. Eventually all times will begin to show convergence. 
This is illustrated by iterations 1 and 3 in Fig. 14. 

There obviously exists a stability for the converged solution in 
the open-loop method. Since the behavior of this stability problem 
was not consistant nothing specific could be found that attributed 
to the problem. The digital to analog output was adjusted to ensure 
no bias voltage existed. The error between the digital and analog 


(1) 


solution for Y was compared and the difference was within the re- 


solution capabilities of the equipment. The differencing of X and 
y was accomplished within the digital computer to nullify any error 
in the analog to digital circuit. The sampling rate at the analog to 
digital interface was varied from 250 to 400 Hz (400 Hz was the maximum 
attainable). In each case the behavior reported for Fig. 16 was 
observed. 

The stability problem appears more pronounced when K.(t) +k, (t)= Kido, 
The neutron kinetics, especially for a fast reactor with Plo x ioe 
are very sensitive to K(t) and if this quantity does not go to zero in 
the limit of steady state then the equations respond accordingly. It 
is a necessary procedure in analog simulation to ensure that with 
K(t) = 0 the reactor kinetics remain in the steady state condition when 
in the compute mode. With equation 3 which requires an integrator with 

4 


a gain of by “= 10 the system requires constant checking to insure a 


balanced condition between N and D so that the steady state solution 
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does not drift in the compute mode. The use of equation (5), which 
removes the integrator with gain Fhe , eliminates this need for constant 
checking, and it was hoped that the stability and convergence problem 
would be reduced. However, no noticable improvement occurred. 

It can be seen fram Fig. 15 that the algorithm developed in this 
model rapidly produces N(t) especially for times less than 1 second. 
However, an error exists for times greater than 1 second between the 
two curves and more iterations should be carried out to reach the 
converged curve. A converged solution is obtained in agreement with 
that obtained in Fig. 14, however, an unusually large number of 
iterations are required when campared with the rate observed in Fig. 
14 in the same time region. Reference 26 reports and demonstrates 
that the rate of convergence is enhanced when starting fram a "good 
initial" guess than fram a "poor" one. It is suspected that the 
stability problem discussed earlier is effecting the rate of conver- 
gence for times greater than 1 sec. More research on this problem 


is needed. 
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VI. CONCLUSIONS AND REMARKS 


The open-loop iteration process as a solution to partial dif- 
ferential equations is a relatively new algorithm for application 
on a hybrid computer. Presently only two reports [Ref. 7 and 24] 
have been located which present investigations into the algorithm. 

It is interesting to note that Ref. 24, which used the open-loop 
method to solve the one dimensional wave equation and heat flow 
equation, addressed itself to results short of the steady state 
solution. Also in Ref. 8, which was a follow up study on Ref. 7, 
the open-loop method was discarded because a large number of iter- 
ations were required near steady state conditions, and a closed-loop 
iteration process was adopted. 

For the open-loop iteration process to be a successful tool in 
solving partial differential equations two problems need further in- 
vestigation, namely, the rate of convergence and the stability problem 
observed in this report. The neutron equations can be reduced to an 
integral operator [Ref. 7] and therefore the convergence is slow 
especially for times near the steady state solution. A large number 
of iterations detract from the advantages of using a hybrid camputer. 
In Ref. 24 the number of iterations in which convergence occurred 
over a given time interval increased substantially as the solution 
approached steady state. 

The problem of stability near steady state solution appears to be 
attributed to the algorithm or the hardware. It is doubtful if the 


model is producing the problem since similar behavior was observed in 
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a very simple system, Fig. 16. Further investigation is needed to 
establish the source of this problem. It is suggested that based 
on the investigation carried out in this report that the problem 
is inherent with the open-loop method. 

One possibility looks encouraging for future studies in solving 
the rate of convergence and the stability problem. The technique 
employed in Ref. 8 utilizes a closed-loop iteration process. The 
method rapidly develops N(t) for later values of time but N(t) for 
intemmediate values of time required more iterations. The algorithm 
in this report develops N(t) for early times but requires additional 
iterations to reach a solution. A cambination of the two may prove 
fruireraan 

In addition the use of the Diode Function Generator represents 
a very flexible tool when used to duplicate a nonlinear function. 
This is especially demonstrated when incorporated in the continuous 
solution to the nonlinear heat flow equation. The procedure used in 
Ref. 7 was to approximate K(T) as a third order polynomial which 
resulted in the need for a large number of multipliers and summer 
amplifiers. Also the DFG can approximate the curve more accurately 
then can the third order polynomial. 

It is also worthwhile to note that the present hybrid simulation 
has the following limitations: 

1) The bonding between the fuel and cladding was neglected. 
This will produce lower temperatures in the fuel than actually en- 
countered. An attempt was made to include the bonding but the 
analog circuit was unstable. 

2) The point model neutron kinetics are used to describe the 


behavior of the entire core. 
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3) Axial flow of heat is neglected. This will result in 
producing higher temperatures. However any error in temperature 
introduced here is considered negligible. 

4) The assumption that an average fuel cell gives all the 
information necessary to obtain the average core temperature may 
lead to unrealistic dynamic behavior [Ref. 6]. However in the 
design of the General Electric LMFBR the fissile fuel is varied 
radially within the core to flatten the power flux. This would 
tend to make the assumption more realistic. 

Presently this model is very camplex and some problems with the 
present model need further investigation before the limitations 


outlined above can be considered. 
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APPENDIX A 


A. MATHEMATICAL PROOFS AND DERIVATIONS 


dee Kaine tae Equations 


The kinetic equations are recast here in the reduced variable 


form. Equation (1) is considered with one precursor concentration, 


namely, 
kale) [ pit) -e] Mle) + > Ce) (4) 
At fia | 
(1) 
Acit) 8 \. \c(Z) b) 
Soe & Mit)-A ( 
G ; 
wheve, = = PAs Re 


Equation (la) is divided by yn(0) and equation (lb) by c(0), 


which are the values of neutron and precursor concentration (#/cm?) 


respectively at steady state. Defining two quantities, N(t) = n(t) 


n{(o) 
and D(t) = — equation (1) becomes 
ANE | Ce 
Sou ees = polls Q 
a7 7g | p(t) BING) + d oe (a) a 


ees eee Noe 
dt 4 (>) } 


d vuldiueeel Si 


At steady state Ty Te =O 


in equation (1), therefore, 


% VCO) eer, CONN. 


Putting this into equation (27) the following 
equations are obtained, 
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UN eae (3) = | 
abet. ~ Le) Blniy + £ Dé) fe) 


) 
- 2 PE net)- nee) + dea) J (28) 
AM [wet d4)] (L) 
“dt 
The quantity OUR = K(t) has the units of dollars. Thus one dollars 
worth of reactivity indicates p (the excess reactivity) has the value 


of 8 . Introducing the term K(t) equation (2) is obtained. 


A NU) / 
aie o | KH) NG) = (4) + DU) (“) 
A DU) , 


TD en as (b) 
where at t = 0, N(0) = D(0) = 
2. Reactivity Feedback 
a. Doppler Effect 
In a LMFBR the Doppler effect contributes the majority of 
the feedback [Ref. 9]. A highly simplified picture [Ref. 25] of the 
Doppler effect is illustrated in Fig. 19. The cross-section curve 
which shows the probability of a nucleus-neutron reaction is shown 
at left. At the right (Case I) is a picture of a neutron approaching 
a nucleus at a low temperature. Since the nucleus is essentially at 
rest neutrons with a lesser or greater velocity than that occurring 
at resonance will escape the reaction. However as the temperature 
increases (Case II) the nucleus is no longer at rest and neutrons of 
lower or higher velocity can still achieve the same relative velocity. 
The effect is to broaden the resonance peak and increase the chances 
of a neutron being absorbed. The absorption of neutrons for other than 


fission processes dominates and a negative reactivity occurs. 
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Figure 19. Illustration of Dopper Effect. 
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A development of an exact mathematical model for the 
Doppler effect has been hampered due to lack of experimental data 
for neutron cross sections at higher energies. It 1s necessary to 
use statistical treatments [Ref. 26] which yield for the Doppler 


effect the expression: 3 
5 =A 
ae ie 
OT 
However recent results by R. Froelich, and others [Ref. 
27] indicate the Doppler effect to be 
300\ 0 
#.-¢( >) 
T 
where C and y are positive constants depending on volume ratio of 
UO., to PO, and sodium density. y is approximately equal to one. 
In the General Electric reactor under design y is equal to one. 


Therefore the Doppler effect is: 


dp = = alos: 
IT = 


where T is the average temperature of the core and Dopp is a position 
constant termed the Doppler coefficient. Starting at an average temn- 
perature of the core, T , the change in reactivity due to Doppler 
O 
effect ( can be expressed as follows: 
Poe) . 


x. 
we 7 diam AT 
ie _ Poser / a Aoep Ay, (1+ = 


7, é 
or, in terms of dollars 


b. Fuel Expansion 
The expansion of the fuel gives rise to a negative feed- 
back. It can be explained somewhat simply by observing that an 
expansion of the fuel reduces the fission nuclei per om” seen by the 


flux of neutrons. In the General Electric design LMFBR it is expressed 


as 
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or, in dollar units 


I3 a 
Kexp = ~ anal aT; 


where a is a positive constant. 
c. Coolant Expansion 

The expansion of the coolant can give rise to either a 
positive or negative feedback. As the coolant expands the number of 
neutrons captured parasitically is reduced thereby increasing reac- 
tivity. Changes occur in neutron leakage (loss of scattering events 
by sodium) when the coolant expands and, hence, reduction of reac- 
tivity occurs. Also less moderation occurs producing a hardening 
of the neutron spectrum, and if fertile materials are present more 
fissions occur producing an increase in reactivity. The size of the 
core, fuel and volume ratio of coolant determine which of the above 
are more pronounced. In the General Electric reactor the feedback 


due to coolant expansion is positive and expressed as follows, 


P, = oo aT. 


OO} 


or, 8in dollar unites, 


Keoo, = Ceoes OT /g 
where C col is a positive constant usually referred to as the sodium 
void coefficient. 
3. Heat Flow Equations 
Fram Section IIA the equation ee flow in the fuel rod 
is, 


ora oT = = 3 | K,(T)y aT | S(2z) N(t) (9) 
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Figure 20 illustrates the various radial segments selected within the 
fuel pin. The temperature between mid-points of each radial segment 
is assumed to vary linearly. The slope at each side of a segment 
interface are equal thus ensuring that the flow of heat across the 
boundary is continuous. Applying equation (8) to the region indi- 
cated in Fig. 20 and integrating over the spacial variable r the 


following equation is obtained,” 


r 


v.44" ir S cas 
t = 
a | K (Tv oT | des y s@)NQ@) 4v 
V.-a¥ —as 
ee ome tc 


For clarity each integral will be handled separately. Assuming K(T) 


is not a constant, the first expression on the right side becomes, 





v +o Vy or 
( Air [ ar | . 
; Sf K, (ry wT dy = La) Nae Pa av 
aed vied 
ee ie 
= KG ais Creston — ) 
, 2")( Ov (29) 
ore 
: 
= gare (7-47) OY ) 
where, K ont = thermal conductivity at the boundary between 


e ® . : 7 
node i+ 1andi. The temperature at the boundary is Tent 


The second integral on the right side reduces to, 


Reference 15 states this technique of reducing a partial dif- 
ferential equation to an ordinary differention equation requires a 
minimum number of nodes. 
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yr 


Vos 
_ 2 a. 
[vscncrae = (2) Nt) a 


Y. a ar 
- 


. (30) 


— 


ela Nor) Vee 


The quantity c,d. in the integral on the left hand side of 


equation (9) is assumed to be a Seneeenie. therefore, cd. and the 
operator Ye 


can be brought outside the integral. The 


integral is then handled as follows: 


cal 


Cre eT ee Cz Sr rane fr | fe) 
[ves x a 
rT oe 


t — 


From Fig. 20 the relation es . 1s, 
Tar ~ (EBV) (et ere 


- (len) weve (i482) 


Therefore equation (31) becomes, 








In most materials this can be considered a constant and in UO 


the variation is not pronounced. Reference 15 treats the conditiof 
where Crd. is not a constant. 
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Vv, + at (32) 





After applying limits and same reduction equation (32) becomes: 


d if iy. ve 2 
= Gd = ie m ae, a eee 
ot 24 4 
(33) 


y 


The operator Y , in equation (33) can be changed to At Since the 
Space variable has been eliminated. Cambining equations (29), (30), 
and (33) and after dividing through by Y,; 4v_ the following equation 


is obtained, 
d, & [t. (4-8, )+ ZT 4+ Te (G- Z,,) 
Ce Ae Lori Ns oe 4 ae dN 
i ayY a coal aY Ti ee 
= eee (1+ \( a ) KU a —') (34) 


+ S(2) Nt) 


The equation for heat flow in the cladding is, 


OT ft oO 
ad Ze t 2 [acne Z| 00 


Ve 


For the cladding which for this model is stainless steel the 
quantities cd. and K, are assumed constant. Following an identical 
procedure as above and utilizing Fig. 21 equation (10) can be reduced 


to, 





cod 4 [a (4-28 )+m , Ta(4r 24 )] 


=k (14 )( Be) Ke Be) 


To differentiate one axial segment from the next, superscript 


(35) 


j is used and equations (34) and (35) become equations (11) and (12) 
used in Section IIB. 

To the clad-coolant interface the conservation of energy is 
applied to develop a differential equation which couples each axial 
segment. As stated earlier the axial flow of heat is neglected and 
therefore no coupling of one axial segment with its neighbor exists 
within the fuel and cladding regions. The coupling occurs through the 
flowing coolant. 

Figure 22 represents a top and side view of an axial segment. 
If the coolant is treated as a heat sink, the temperature drops at 
the clad-coolant interface to a value assumed constant (Fig. 21) over 
the region which receives heat from the segment in question. Ar col 
represents the extent in the radially direction of this region. The 
radial gradient of temperature in coolant is neglected. 

The heat (éq) received by the coolant is 


bg a nil Ye (a2) Re) Ceou) ne able (t,2) 
(36) 


oe cl ie (1,2) = AL awe: oe dt 
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Figure 22. Diagram of Axial Segment. 
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therefore equation (36) becomes 


ee dec. | SE 42 r cE at] 








Zila 2 Gos, (33),) 
or, 
dz 
j Jd 4 | =. = aie dt ai 
SLOOP TNS VIA 2 ee) (38) 
; eee Jt 
azo 7 
aaae Mass flow rate of coolant = G(t) 
dg | 
—————— = Heat flux across clad-coolant interface = H 
2h Croker G 


Both numerator and denominator of the left side of equation (39) 
are now multiplied by 27x. in order to cast it in a form to take ad- 
vantage of certain design criteria. In a reactor the volume fraction 
of coolant is a known design criteria and this knowledge allows the 
evaluation of 2nx, A which is approximately the cross section 
of coolant. In most LMFBR designs [Ref. 28] considered to date, includ- 
ing the General Electric model used herein, the volume fraction of 
coolant <50% implying the cross section of coolant equals the combined 


cross section of fuel and clad. Therefore, 


Zien 
rer eae ratio of circumference of clad to cross sectional 
TY. OX eo, area of coolant 


DT Ns. 
age 





Equation (38) may now be written as, 





Gt) Wm , oT 2H 
ees Oe ot a aoe crn V, (39) 


dS 


Again an ordinary differential equation is sought in order to 
use the analog computer. To reduce equation (39) to an ordinary dif- 
ferential equation the finite differencing technique [Ref. 29] is 


employed by discretizing the space variable. The operator Yo is 

Tatas 
at 

T. across the Az segment, thus, the average coolant temperature (T) 


reduced to Also, a linear relationship is assumed for 





for a segment can be expressed as 


=e eT + Tae 
. 7a 
j= 1, 2. °10 (10 being the number of segments) 


Ten axial segments were found adequate to describe the coolant ten- 
perature profile in agreement with Ref. 7. ie represents the outlet 
temperature from the previous segment and the outlet temperature of 
the segment in question. 
With these substitutions equation (40) becames, 

AT me Gt) aig mete, wii 

_—_— 5 ee ee 
but, 1 = ge ca 


therefore, 





—— tf 





Ape Zea ( 
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ry 
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Coot le, Coby eG (14) 
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APPENDIX B 


A. SCALING OF DIFFERENTIAL EQUATIONS 


The scaling of the differential equations (15) through (21) 


will be demonstrated by using equation (15) 


a [7 + .26337% |= Be eyes LK, 3 (T'-T')] +. 4100 $(2) N(t) (15) 


The various quantities (T, 1 Tor Ke /2 etc.) are represented on the 
analog computer as voltages whereas in equation (15) they represent 
units (CK, watts/an”, etc.), therefore a coefficient called a scaling 
factor is needed which transforms the various units into voltages. 


Thus the following expression for each variable in equation (1) is 
defined. 


R= ee 
The ke he 
Kez Re ae 
Ni)= kb, NG) 
S@)= Rk, $(2) 
where Th T., Ke /2 N(t) and §(z) are in volts. Placing these 


quantities into equation (15) and dividing through by kp gives 


dts sasw Je 3202 [ee %, GT] 
(40) 
+ .4] ( ) s(2) NU) 


The value of the scaling factors Kye Ke Ke and ky are determined by 


knowing the maximum voltages the variables will assume during the 
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solution. These are sametimes not known and often the problem must 
be put on the computer and tested before arriving at scaling factors. 


In this case it was found convenient to use the following values for 


the scaling factors, 
fn = ioe /aan 


. _ Watts 
K ~ Cym-°C- Volt 


Watts 
= 00 ——— 
Rs | Cyn3= Volt 


Ry = * Svoit 
With these scaling factors equation (40) becomes 


4 [A +. 26337, ile 32,82] 10 Kea) iy 7 7.') +205 5(2)N) (41) 


Similar operation on the remaining equations (16) through (21) produce 
the following equations 


4. ie Sp EOE om 1295 7, | = 23. iL1o Kiar (Ty- T')| 


= (42) 
Sess [10K 9, (F-T)] +. 2165 £@2) N Ct) 


; a _ oe 
+, ire +. 1889 Ta +. 4447s | = 20.8 [io kK, 4, (T- 13 )] 
(43) 


— /3.88| 10 Kis. (1 -Ty ) + 21657 $@)NQ) 


a eae ek: —; : 
AW surT ents le m0 FER 
(44) 
— 14.85 fio Ky 4, (Ti -T) J 4 Lies 3 (2) NCt) 
d f=; = a _ 
Ge LTE + «129074 +4362 7 [= 46.7 [ro Ky, (-)] 
(45) 
—1875 [10 Ky, (Te -T)) +.264 Sa) NU) 
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4. (Ti+ /.ses T+ east |= 174.1 [R- T ] 

~ 69.8 [io Ry, (T- Ts')] 
on 3 _ =; ae! ose 
Ae -~t5es te + 680, + 5370aee 


(Note: ie refers to average coolant temperature in volts) 
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APPENDIX C 


A. PHYSICAL PARAMETERS (7) 


~ 


Specific heat-c (watt-sec/gm-—C 
Fuel - 0.3345 


Cladding- 0.502 
Coolant - 0.21 

Thermal conductivity-K (T) (watts/am—"C 
Fuel - see Fig. 5 


Cladding - 0.2075 


Density-p (gm/em?) 
Fuel - 9.2 


Cladding - 8.0 
Cogilanic Unis 
Mass flow rate of coolant-G(t) (gm/an*-sec 
500 
Heat transfer coefficient between clad and coolant-h (watts/am?—°C 
13302 
Dimensions (cm) 


Ar = 0.05 


Ar = 0.05 
c 

Az = 10 

L = 100 


Prompt neutron lifetime-2 (sec) 
4.8 x 107? 


Values taken from Ref. 9 and personal correspondence with Mr. C. K. 
Sanathanan. 
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Decay constant for six precursor neeiali(secns) 


i ai 
1 0.0130 
2 0.0314 
5 OF136 
4 0.340 
5 Bl oe/9, 
6 3.500 


Delayed neutron fraction from jth group~8i 


i gi 
1 -0000759 
2 ~000626 
S .000564 
4 -0010700 
5 -000489 
6 .000163 


G 
p= 2 Bi = .002988 


Decay constant for one precursor psa oleiisecn>) 
je APs seas 
Ajopp = +005 mm 
be 
B = eaeeenin 
exp : 
C oey = 2x 10 
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